REPOIIT  DOCUMENTATION  ^AGE 


i«.  RCPORT  sccumrv  cuts 


2».  S6CUWTY  CLASSl»<ATiO*«  AuTHO«»rv 


2b.  OCCLAS$i^lCAriOM/OOW*iCJ)iAO*MC  X>«iOvU 


4.  PERFORMINC  OROANUATIO**  REPORT  buMRiPil] 

ME^  Tectmical  MenDrandun  401 


6a.  NAME  OF  PERFORMINC  ORCANtZAnON 

Marine  Physical  Labcoratory 


i  mCi'u.'oiRiiw;  iON 


«•  04F<I 

OP  ^IlgMRWI 


8a.  NAME  OF  FUNDING /SPONSORINC 
ORGANIZATION 

Office  of  Naval  Research 


8c  ADDRESS  (Cty,  Start,  aod  ZfPCodtJ 
Department  of  the  Navy 
800  North  Quincy  Street 
Arlington,  VA  22217-5000 


1 1 .  TITLE  (Indudt  Setunty  Q^aihcstioni 

PARABOLIC  BQUATIGN  MODEL 


12.  PERSONAL  AUTHOR(S) 

Jean-Mane  Tran 


13a.  TYPE  OF  REPORT 

Sumnary  _ 


16.  SUPPLEMENTARY  NOTATION 


to  04<<f 

0*  a#0*<a»itj 

CNR 


ASiSaRtl^  i<.t»  iat**  J'lfCami 

8CO  ?4esrth  Ouine>* 
Arlingtan.  VA  -5-000 


*  taO(  PWt'RtiWlttW*  ICA*  ICIK 

}iOOOU-«7<-QJJ7 


0  VO-uKE  O*  < 


pRCCAaw  I  toC  4  <  * 
I.4V4NT  NO  I  NO 


(VOP*  UNit 
A<Cll.ViON  NO 


13b  TIME  COVERED 

FROM  TO 

34 

1 

t*0*'  jhetr  ManaWi  Der/  1 

1388 _ _J 

17. 

COSATI  CODES  1 

FIELD 

GROUP 

SUB-GROUP  1 

19.  ABSTRACT 


'8  SuBjEC^  terms  (Contanvt  o**  -T  a«t  tN  b*ot* 

parabolic  propagation  nodel,  acoustic  rtxSellino, 
SACLANT  PPifBQ  prograro 


This  report  deals  wia  a  pibolk  pwpNtMW  awM  at*  aaplHaBMM  «  d* 
Maiiiie  Physical  LabocMOty.  This  panbabc  tRaMa  (K)  profraai  a  s  FORTR.A>'  tt 
venioa  of  (he  SACLANT  PAREQ  pfosw  (U- 

ThepMaboUcweihodiwtwaiaaior  dnelopwiwraidNfaMofatqaaat  aaodHiiv 
[2],  The  sdvaniace  of  die  aiediod  if  dia  simpla  way  laaft  depcadeaGC  t  rinilliil  tM- 
fiacdon  eflecB  are  included,  but  netiber  irrcrbcrwoa  aor  backsoMcnaf  m  matmmi 
for  (31.  The  pmbolic  equMion  lequiias  mall  Mslca  at  boat  m  Huruoaui 

and  cannot  handle  steep  sound  vetociiy  padicaa  or  Mgb  ftaqaeacy  pnpa^sMi  becaaK 


of  computer  speed  considetanona.  The  nauung  tiwa  is  pwpttwml  lo  W/  ”  —  R_ 

I  •* 

where  N  is  the  water  plus  boaom  heiiht.  /  the  frequency,  c  die  sonad  wetacny.  t  the 
depth  and  maximum  ran|e  (4). 

In  the  following  peges.  the  pmbolic  equations  m  denved.  die  ataocuaed  anmen- 
cal  algorithms  are  presented,  the  enviiamnem  end  die  progiam  iaputs  are  described. 
Finally  examples  are  given. 

The  goal  of  this  report  is  to  provide  the  potemiel  user  the  necctmy  mformauan  to 
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ABSTRACT 


This  report  deals  with  a  parabolic  propagation  model  now  implemented  at  the 
Marine  Physical  Laboratory.  This  parabolic  equation  (PE)  program  is  a  FORTRAN  77 
version  of  the  SACLANTPAREQprogram.[l].  r 


The  parabolic  methods  were  a  major  development  in  the  field  of  acoustic  modeling. 
[2].  The  advantage  of  the  method  is  the  simple  way  range  dependence  is  handled.  Dif¬ 
fraction  effects  are  included,  but  neither  reverberation  nor  backscattering  are  accounted 
rorT3I:  The  parabolic  equation  requires  small  angles  of  propagation  firom  the  horizontal 
and  cannot  handle  steep  sound  velocity  gradients  or  high  frequency  proparation  because 

of  computer  speed  considerations.'  The  running  time  is  proportional  to 


dz 

where  //  is  the  water  plus  bottom  height,  /  the  frequency,  c  the  sound  velocity,  z  the 
depth  and  the  maximum  range  [4]. 


In  the  following  pages,  \he  parabolic  equations  are  derived,  the  associated  numeri¬ 
cal  algorithms  are  presented,  the  environment  and  the  program  inputs  are  described. 
Finally  examples  are  given. 


The  goal  of  this  report  is  to  provide 
understand  and  use  the  PE  program,  f  ^ 
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1.  Parabolic  Approximation 

The  starting  point  of  acoustic  models  is  the  homogeneous  Helmholtz  equa¬ 
tion  or  reduced  wave  equation  which  is  an  elliptic  partial  differential  equation  : 

7^p+k^p=0  (1.1) 

V^p+k^n^p=0  (1.2) 

with  k^  =  kgn^  where  k^  is  a  reference  wave  number  and  n  can  be  regarded  as  a 
refraction  index.  In  cylindrical  coordinates,  neglecting  the  azimuthal  depen¬ 
dence,  the  Helmholtz  equation  is  : 

Intrinsic  cylindrical  spreading  is  accounted  for  by  a  factor  in  the  pressure, 
leading  to  the  classic  change  of  variable  p  =  r"’^(r  ,z). 

Then,  the  equation  becomes  : 

-^+-p-K*o^n2+— )«1>  =  0  (1.4) 

Essentially,  the  far  field  zone  is  of  interest ;  therefore,  k„r»l  is  assumed  with  an 
index  of  lefiraction  of  the  order  of  unity.  The  equation  is  simplified 


where  the  operator  Q  is 


(■^+^/C>=o 


^2  -2.  1  o 
^  ■  k^  hz^ 


Parabolic  approximations  are  obtained  by  factorization  of  the  differential 
operator  Q  into  the  product  of  two  commuting  operators 


i 
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=  Hj;+k,Q)(ijp*k,Q)  ( 1 .8) 

It  requires  that  and  Q  commute.  This  condition  is  considered  saiisAed  if  the 

index  of  refraction  slowly  varies  in  range  so  that  the  gradient  ^  is  negligible. 
The  differential  equation  then  can  be  separated  into  two  equations  : 


{if:+k.Q)O  =  0 


(1  10) 


The  general  solution  is  the  sum  of  their  two  solutions,  the  first  one  represenu  the 
outgoing  wave  and  the  second  one  the  incoming  wave.  The  underlying  nature  of 
the  parabolic  approximation  is  to  consider  that  the  solution  propagates  primanly 
in  the  outward  direction  with  very  little  reflection,  backscattering  and  leverbera- 


■i 
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2.  Parabolk  Equations 

The  differences  in  formulations  of  the  parabolic  approximations  leading  to 
din’ereni  parabolic  equations  reside  in  the  way  the  Q  operator  is  evaluated.  The 
pseudo  differential  operator  Q  is  rewritten  as 


(2.1) 

ko 

and  linearized  in  different  ways  by  assuming  either  (n^-l)<l>  or  -^4^  is  small 

Kq  oz 

compared  to 

Taking  is  equivalent  to  considering  small  gradients  of  the  sound 

velocity.  Tapper!  (6]  shows  that  the  term  \  is  related  to  the  mean  square 

ie  OZ 

angle  of  propagation  with  respect  to  the  horizontal.  If  <l>  =  exp  ik^  (r  co‘'>0±2  sin0),  the 
norm  of  -4^^  ‘s  where  6  is  the  angle  of  propagation  with  respect  to  the 

horizontal.  Therefore  assuming  it  is  small  is  equivalent  to  considering  small  pro¬ 
pagation  angles  with  respect  to  the  horizontal. 

Approximations  of  the  Q  operator  lead  to  two  classes  of  models,  one  of 
iterative  methods  based  on  finite  difference  schemes,  the  other  based  on  the  Split 
Step  Algorithm  (6,7]  using  Fast  Fourier  Transforms  (FFTs).  Only  the  last  class 
will  be  descibed  here  in  relation  to  MPL’s  PE  program. 

The  standard  small  angle  parabolic  approximation  is  obtained  by  assuming 
that  both  the  index  of  refraction  term  and  the  propagation  angle  term  are  small 
and  by  using  a  first  order  Taylor  series  [5,6]. 


^  I,  1  ^ 

2  ^2k/dz^ 
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The  corresponding  equation  is  obtained  by  the  classic  transformation  «l>  =  : 

(‘‘*‘'^  =  0  (2.3) 

^+2iit<,-^+Jk,2(n2-l)'P  =  0  (2.5) 

Equation  (2.5)  is  the  standard  parabolic  equation,  also  called  Tappen  and  Hardin 
parabolic  equation. 


Thompson  and  Chapman  [7]  introduced  an  improved  parabolic  approxima¬ 
tion  by  using  another  development  of  Q  proposed  by  Feit  and  Reck  |8) ; 


They  showed  that  this  approximation  can  stand  much  wider  angles.  While  the 
standard  parabolic  equation  is  valid  to  about  ±20*,  this  modified  parabolic  equa¬ 
tion  is  valid  up  to  ±40®  [4].  Using  an  additional  pseudo-operator  identity 


(l-HiV'^+l 


the  equation  becomes 


±JL 

1+ - ; — rj - +n-\  «&  =  0 


With  the  previous  envelope  transformation  «I>  =  one  gets 

!  3^ 


(2.10) 


X  Numerical  Iropkinciitalioo 


These  {MuaboUc  equations  are  extensively  used  because  of  their  efficient 
numerical  implementation  based  on  FFTs.  Jensen  and  Krol  [2,9]  elegantly  derive 
the  Split  Step  algorithm  for  the  Tappen  and  Hardin  parabolic  equation.  They  do 
us  Founer  transform  and  assume  an  almost  constant  index  of  refraction.  By 
inverse  Fourier  Transfwm.  they  get 


2  f-l 


e  F(4'(r„2)) 


(3.1) 


where  f  is  the  one  dimensional  Fourier  transform  with  respect  to  z,  is  a  refer¬ 
ence  range,  and «  is  the  Fourier  variable. 

The  above  algorithm  was  primarily  derived  by  Tappert  and  Hardin  [5,6] 
using  operator  represenutions.  The  parabolic  differential  equation  is 


~=iA'V+iB'¥  (3.2) 

dr 

It  13* 

with  j)  1)  and  B{2)  =  — - r  .  Note  that  at  this  point,  attenuation  can 

2  2kf  dz 

'y 

be  introduced  by  inserting  a  complex  term  in  a  :  a  =  -|-tn‘-i+ia).  The  solution 
is  then  assumed  to  be 


4'(r,+Ar,z) 


exp 


i  J  A{rj)dr+i  j  B{z)dr 


'?('•<,  .2) 


(3.3) 


[  ['■  '•  JJ 

by  analogy  with  differential  equations  in  the  real  space  with  constant  coefficients 
A,B.  Assuming  that  A  varies  slowly  with  respect  to  r,  on  splits  the  operator 
'f'(r,+Ar,2)  =  cxp(iAri4)cxp(iArfi)4'(r„,z)  .  The  term  exp{iArB)'¥{rg^)  is  evaluated 


via  its  Fourier  transform 


exp{iArB)'V(r,;z) 


(3.4) 


4.  Interpretation  of  the  Split  Step  Algorithm 


Two  physical  phenomena  are  involved  :  the  propagaiton  in  a  homogeneous 
environment  and  the  influence  of  changes  in  the  medium  on  the  pfopagauon 
The  Split  Step  Algorithm  multiplies  the  angular  spectrum  by  a  phase  shift  to 
account  for  propagation  over  a  range  step  jv  The  steeper  the  angle,  the  greaier  is 
the  phase  shift.  Then  it  muluplies  the  inverse  Founer  transform  of  this  nca  range 
step  spectrum  by  the  index  of  refraction  term  which  accounts  for  changes  m  the 
medium. 

Other  algorithms  for  the  standard  parabolic  equation  were  introduced  |h.)0] 
such  as 


y  (r,  j )  ■  eapi 


'■'1 


cap 


2*. 


f  ittf  t 


It  allows  for  a  better  coupling  between  the  two  physical  phenomena  dcscnhrd 
above.  An  error  analysis  ( 10)  shows  that  the  tninal  algorithm  (3  )  i  has  a  second 
order  accuracy  in  t  while  the  last  one  (4  I)  has  a  third  ortder  accuracy  in  * 

The  Split  Step  Algonthm  is  a  marching  solution  in  range,  ihcrefore  range 
dependence  of  the  environmental  data  f sound  velocity  profiles.  Kaihytnctry »  is 
easily  introduced.  Boundary  conditions  arc  tmpiicii  in  the  numerKal  solution  A 
pressure -re lease  surface  at  the  air/watcr  interface  is  assumed  .  the  field  is  made 
anti-symmetrical  in  depth  about  *  =  o  by  using  real  sine  transforms  ( 10| 


The  Founer  tlteory  as-sumes  L '  integrable  functions,  the  field  must  vanish 
below  the  maximum  depth  of  the  transform.  A  pressure  release  bottom  must 
correspond  to  the  maximum  transform  sample  to  avoid  aliasing  while  doing 


An  unphyturai  butUMn  Uycf  u  uiauduLcd  lo  Aitmujtu;  ihe  ticUl  in  iKc  buM 
torn  tlO).  Thi.%  iechnii)uc  ituu  4t  Umt  (be  liununjint  contnbuiion 

comn  from  U>w  onkr  moiJct  (Km  cont^pund  lo  ioi*  »n$\c  The  higher 

(he  mode,  (he  gtt^uei  it  (he  frmeuig  angle,  and  (hr  more  prnrtradon  iniu  the  hoi 


>?v.v,v 


5.  Dcscriptkm  of  MPL’i  PE  Procrom 


The  atgohihtm  used  are  the  Split  Step  Algorithms  (3  1)  and  (3.5)  which 


have  a  second  order  accuracy  in  r .  The  inioal  transform  size  is  determined  by  the 


input  transform  component  H,  the  transform  has  a  total  of  "f  putnis.  s  must  be 


smaller  than  12.  The  maximum  depth  of  the  transform  is  the  maximum  depth 


of  the  bottom  extended  by  a  factor  33  ^  to  account  for  the  unphysical 


attenuating  layer. 


If  the  input  transfonn  size  is  not  given,  it  is  determined  by  the  program  using 


N  =  U\og^kZ^n.)  (5.1) 

where  k  is  the  average  wave  number,  and  is  the  maximum  depth  transform. 


The  program  begins  by  creating  the  starting  field  and,  in  a  range  loop,  reads 


the  environmental  data,  creates  the  refraction  index  table,  e'^^,  and  the  second 


derivative  table,  .  If  the  range  step  is  not  fixed,  the  program  computes  a  new 


range  step  at  each  range  :  hr  =  - - -  where  X  is  the  wavelength  and 


0  =  tan" 


T  iptxnsi  *  ^  the  Fourier  variable,  F  is  the  Fourier  transform  and 

Kp  1 1*  \  *  )  I 


isthct,*  nonn(i.e.  i/’i  »  \FF'di). 

If  ihe  nov  ranse  step  has  a  relative  difference  with  respect  to  the  old  range 
step  of  at  least  25%  then  this  new  range  step  is  selected  and  a  new  table  of  second 
derivatives  is  computed.  Otherwise,  the  old  range  step  and  the  second 
derivative  table  are  kept.  The  program  does  Founer  sine  transforms  and  multi¬ 
plies  the  transformed  field  by  the  tabte  of  second  denvauves  : 


e*P  -* - ; - ■57“  {oamdardFE) 

•  I 


cap  -I- 


2a, 


FE) 


The  inverse  sine  Founer  transform  is  performed  and  the  field  is  then  multiplied 


e»p|i — #■*“'**  (jM^arrf  FE) 

f  ar«, 

e*p  I— - 1)  r^"  FE ) 


a  is  the  attenuation  in  the  water  or  in  the  bottom.  The  program  advances  the 
range  to  r  > 

After  bookkeeping  on  the  range  dependence  of  the  environment  is  com¬ 
pleted.  it  checks  for  aliasing  in  the  transform.  Aliasing  occurs  when  the  contribu¬ 
tion  of  steep  angles  in  the  angular  spectrum  is  significant.  The  parabolic  approxi¬ 
mation  is  based  on  the  ’’small  angle”  hypothesis.  The  angular  spectrum  is 
divided  into  four  equal  bins,  their  energy.  £«  from  low  to  high 

angles  are  computed.  An  aliasing  factor  is  defined  as  tOlog,o£i  with 

£4 

^1  =  1: — e — e — -14  dS  </?],  then  the  energy  is  spread  throughout  a  large 

C  l+£  2’fC34-C4 


Ljl*'*.i*a^**i 


ia,* ■*4' ■' 


angular  interval  and  aliasing  is  severe  so  the  run  is  terminated.  If 
-20 dB  <Ri<-l4dB,  then  a  warning  is  issued  and  after  five  warnings  the  run  is 

£  ^  £  j 

terminated.  Two  oversampling  factors  are  defined,  /?2=  -= — ^  and  Rj=—.  If 

C  i+C  2  £2 

/?2<-70de  andl?3<-60dB,  the  field  is  oversampled  and  the  transform  size  is 
reduced  by  a  factor  of  2. 


! 


6.  Phase  Error  Corrections 


The  parabolic  approximation  and  the  Split  Step  Algorithm  are  applicable  if 
the  index  of  refraction  locally  varies  only  slightly  and  if  propagation  is  limited  in 
a  narrow  aperture  about  the  outward  range. 

Comparison  between  the  parabolic  equation  and  elliptic  equation  normal 
mode  phase  velocities  in  range  independent  waveguides  [11,12]  showed  that  the 
parabolic  approximation  introduces  phase  errors. 

In  the  case  of  a  single  normal  mode  describing  the  propagation.the  phase 
error  is  canceled  if  the  reference  wave  number,  ,  is  chosen  to  be  the  mode  wave 
number  [11].  When  several  modes  are  propagating,  there  is  no  obvious  a  priori 
choice  for  [13]  even  if  it  can  be  intuitively  understood  as  a  weighted  sum  of 
the  excited  mode  wave  numbers  [12].  In  practice,  the  choice  of  k^  is  left  to  the 
model  user.  Two  different  methods  implemented  on  MPL’s  PE  program  try  to 
correct  these  phase  errors. 

The  first  method  by  Brock  et  al  [14],  also  called  the  CMOD  correction, 
modifies  the  environment,  i.e.  the  profile  of  the  refraction  index  to  minimize  the 
phase  errors.  The  idea  is  to  decompose  the  sound  field  calculated  by  the  para¬ 
bolic  equation  into  propagating  normal  modes  and  to  compare  these  modes  with 
those  obtained  by  a  WKB  approximation  to  the  elliptic  wave  equation.  Requir¬ 
ing  that  the  depths  of  the  modal  turning  points  (where  the  mode  amplitudes  are 
the  largest)  are  the  same  in  the  parabolic  and  elliptic  cases,  that  the  phase  veloci¬ 
ties  at  the  turning  points  are  the  same  and  assuming  that  there  is  no  isovelocity 
region,  a  mapping  can  be  derived  : 

(n,z)^(Vn^-l,2Vjr)  (6.1) 


The  CMOD  correction  appears  to  offer  improvements  when  the  source  and 
the  receiver  arc  near  the  surface,  the  bottom,  or  in  regions  near  individual  mode 
turning  points. 

Pierce  introduced  a  better  approach  to  determine  in  a  more  natural  way  the 
reference  wave  number  by  using  the  Rayleigh  Principle  :  "if  the  propagation  is 
progressive  such  that  the  energy  is  being  transported  on  the  average  in  one  diiec- 
tion  without  propagation  in  the  backward  direction  then  the  average  kinetic 
energy  is  equal  to  the  average  potential  enCTgy"  [13]. 

The  kinetic  energy  in  the  volume  element,  V^,  is  ypi  lul^+  Ivl^)  V,,  where 

u.v  are  the  particle  velocities  in  the  r.2  directions  respectively.  The  potential 
energy  dV  where  P  is  the  pressure  and  dV  is  the  change  of  volume  of  the 
volume  element  with  varying  pressure.  If  the  volume  element  dimensions  are 
small  compared  to  the  wavelength,  the  plane  wave  approximation  can  be  used 
p  V  dP 

and  ^  =  - r)  and  therefore,  dV  =  — Then,  the  potential  energy 

pc*  pc 

becomes  ^  [15]. 

2pc* 

Using  a  volume  element  integrated  over  z,  the  Rayleigh  Principle  becomes 

/-ip(  lu|2+  (v|2)*  (6.2) 

This  principle  can  be  applied  in  underwater  acoustics  where  the  field  is  a 
superposition  of  modes  [13].  The  equation  becomes 

f(— )y  ‘  \<b\^dz  -fp-'  lO,  I 

_ : _ 

Jp->l<&l^dz 


since 


(6.3) 


This  wave  number  is  range  dependent  and  can  be  proved  to  be  range 
independent  in  range  independent  media.  The  method  also  allows  discontinuities 
of  density.  The  natural  choice  for  k,  as  given  by  Pierce  consistently  yielded  more 
accurate  results  than  when  other  schemes  were  used  [13] .  According  to  [1],  it  is 
strongly  recommended,  to  ensure  maximum  accuracy  of  the  PE  calculations,  to 
use  this  phase  velocity  correction  together  with  the  modified  parabolic  equation 
(Thomson  and  Chapman).  MPL’s  PE  program  computes  every  ten  range  steps. 
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7.  Source  Field 

Initial  fields  must  be  computed  to  start  the  marching  solution  in  range.  The 
wave  equation  to  solve  is 

where  is  the  pressure  at  one  meter  [6].  po  is  set  to  unity. 

The  starting  field  can  be  either  created  by  source  routines  available  in  the 
program  or  entered  as  an  input  in  a  file.  The  three  available  sources  are  the 
Gaussian  source  [6],  the  Greene  source  [16]  and  the  Thomson  source  [7]. 

The  Gaussian  source  is  an  asymptotic  approximation  to  spherical  spreading 
in  the  far  field.  The  standard  parabolic  equation  is  simplified  in  the  near  field  by 
neglecting  the  refraction  term ; 

with  p=r~^e*’''V.  Near  the  point  source  at  depth  z,,  the  field  is  spherically  spread- 

ik.R  _ 

tttg  P-—^  with  R  =  Vr^+(z-r,)^  To  do  an  asymptotic  matching  between  the  two 
R 

expression  for  p,  one  lets  r-x**  in  the  spherical  spreading  solution  and,  after  a 


first  order  expansion 


1  1  (z-z,? 

p  =  -e  ^  1-- - — 

r  2 


Taking  the  solution  of  the  simplified  parabolic  equation,  one  lets  r-»0,  and  by 


asymptotic  connection  one  gets 


'¥(0.2)  =  ^e  ^  (7.4) 

This  solution  given  by  Tappert  is  simple  and  avoids  spurious  sidelobes  [5,6]. 
This  source  has  been  modified  to  allow  some  tilt  with  respect  to  the  horizontal. 
The  half  width  01  of  the  source  (half  angle  corresponding  to  a  3  48  loss  with 


respect  with  the  center  of  the  source)  can  be  selected  too.  Taking  0i  =  0 
corresponds  to  the  original  source  or  the  standard  Gaussian  source,  that  is  to  a 
half  width  of  35** . 


The  Greene  source  is  documented  in  [16].  It  has  been  designed  for  a  higher 
angle  parabolic  equation  (HAPE).  It  is  a  wide  angle  source  with  a  -3  dfi  half 
width  of  80® . 

The  Thomson  source  is  an  alternate  starting  field  that  uses  a  Kaiser  Bessel 
window  to  design  a  source  that  insonifies  a  limited  aperture  in  the  k,  space.  The 
parameters  describing  this  source  are  the  depth  of  the  source  r, ,  the  half  width  of 
the  source  aperture  e,,  the  tilt  of  the  source  02  and  a  the  level  in  dB  between  the 
stopband  and  the  passband  in  the  k,  space.  The  field  is 


'P(0,r)  =  (-r)’^ 


_  i_sin[2rt-^(z-z,)tan0,]-^e 

/  n(z-2,)  C,  laix) 


(7.5) 


with  z:  =  0.1102  (a-21)  and  a>  50.  The  Thomson  source  as  implemented  in  PE  uses 
a =60  so  z:=  4.2978.  I.  is  the  modified  Bessel  function  of  O'*  order. 


y  =  (i-[  ^  ^  1^)''^  where  is  the  maximum  transform  depth  and  c»  is  the 
reference  sound  speed. 


The  PE  program  can  also  be  started  by  a  source  field  stored  in  a  file.  The 
file  must  have  2^  pairs  of  the  real  and  imaginary  parts  of  the  pressure.  Each  pair 
corresponds  to  the  pressure  field  at  a  point  of  the  mesh  which  has  a  depth  incre- 

ment  equal  to  .where  Z„„  is  the  maximum  depth  of  the  transform.  This 

field  should  be  normalized  such  that  the  field  corresponds  to  a  unit  pressure 
source  at  one  meter,  and  be  a  solution  of  (7.1).  The  Gaussian  source  and  the 
Thomson  source  have  been  properly  normalized. 


8.  Description  and  Treatment  of  the  Environment 

The  treatment  of  the  environment  is  described  in  the  documentation  of  the 
SACLANT  PE  program  [1].  The  bottom  is  a  two  layer  fluid  bottom  :  a  sediment 
layer  characterized  by  a  height,  an  arbitrary  sound  velocity  profile,  a  density,  and 
a  compressional  wave  attenuation,  and  a  subbottom  layer  characterized  by  a 
sound  velocity,  a  density  and  a  compressional  wave  attenuation. 


C^(z) 


SEDIMENT 


WATER 


Cj^(z) 


P, 

It* 


SUBBOfTPOM 


The  change  of  density  between  the  various  layers  (water,  sediments,  subbot¬ 
tom)  is  taken  care  of  by  the  method  given  by  Tappert  [6].  The  correct  wave  equa- 


uon  IS 


pV(lvp)+-^/7=0  (8.1) 

P 

where  p  is  the  density.  By  doing  the  change  of  variable  ?  =  the  equation 


becomes 


V^q+kon^q  =0 


where  k„=—  and  n  is  an  effective  index  of  refraction 


K* 


(8.3) 


The  discontinuities  of  the  density  are  not  allowed  and  the  program  smears  out  its 
jumps.  If  the  density  jumps  from  pi  to  P2,  the  smoothed  density  is 


,  ^  (P1+P2)  (P2-P1) _ 

p<2) - ) 


and  L  is  taken  such  that  L  =  — . 


f 

(8.4)  ! 


The  subbottom  has  a  height  twice  the  wavelength.  If  the  phase  velocity  of 
the  subbottom  is  set  to  zero,  this  height  is  also  zero  and  there  is  no  subbottom. 


Bottom  loss  is  introduced  by  adding  a  small  imaginary  part  to  the  wave 
number  [2],  k  ^kg+ia  with  where  kg  =  This  attenuation  is  expressed  as 

part  of  the  refiraction  index  assuming  that  a«— .  The  ima- 

c  C  (D  c 

ginary  part  of  the  refraction  index  is  then  expressed  in  term  of  the  real  part 
n^=^—ng.  The  actual  input  for  the  attenuation  is  p  m  4S/X  where  X.  is  the 


wavelength.  Then  a  =  p 


(0 


-  and  n?  =  - Hg, 

2rtc(201og«)  '  27.287527  * 


Volume  attenuation  in  the  water  is  introduced  by  multiplication  of  the  field 
by  a  term  [1]  but  this  term  is  fairly  low  at  the  fi-equencies  of  interest  for  PE, 
that  is  below  200  Hz . 


Since  the  environment  can  vary  in  range,  one  must  provide  the  program 
with  sets  of  the  environmental  data  at  specified  ranges.  A  sector  is  defined  by 
these  '»ges  associated  with  environmental  data.  In  each  sector,  the  sound  velo¬ 
city  profile  is  range  independent  while  the  bathymetry  is  linearly  interpolated  in 
range  between  the  two  ends  of  the  sector  (water  depth,  sediment  depth,  subbot¬ 
tom  depth).  The  sound  velocity  profile  is  interpolated  linearly  with  depth.  How¬ 
ever  the  sound  velocity  profile  will  jump  from  one  sector  to  another  because  it  is 
not  interpolated  in  range. 
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9.  Program  Inputs 

The  program  has  a  readable  input  file.  Data  can  be  divided  into  two  groups. 
One  has  the  general  information  such  that  frequency,  source  depth,  filenames  for 
output  results.  The  other  contains  environmental  information  that  must  be 
repeated  for  each  sector.  A  list  of  the  inputs  follows.  All  of  the  following  infor¬ 
mation  must  be  included  or  lines  left  blank. 

Al.  Title  of  the  input  file  or  run 

A2.  five  output  filenames  of  seven  characters  each  where  the  output  will  be 
dumped  (sector  information  file,  transmission  loss  versus  range  at  selected 
receiver  depths,  transmission  loss  versus  depth  at  selected  ranges,  SIO  data  file 
for  subsequent  use  for  contour  plot  and  SIO  outfilename  for  the  pressure  field).  If 
blank,  no  file  is  generated. 

A3.  Transform  size  (integer  smaller  or  equal  to  12),  0  meaning  that  the  pro¬ 
gram  itself  will  compute  the  transform  size  based  on  a  sufficient  sampling  for 
the  Gaussian  source. 

A4.  Reference  sound  speed  (allowing  the  calculation  of  /fc„).  A  negative 
sound  speed  means  that  the  Pierce  phase  correction  will  be  used.  A  nul  sound 
velocity  means  that  the  CMOD  correction  will  be  used.  Otherwise  the  given 
sound  velocity  will  be  used  as  reference  sound  velocity  in  the  calculation  of  the 
index  of  refraction,  and  its  value  must  have  a  physical  value  (e.g.  1500  m/s). 

A5.  Output  range  step  size  (in  m  ).  It  must  be  greater  than  the  range  step. 

A6.  Frequency  (in  //z)  since  at  high  frequency  the  range  step  is  small,  the 
frequency  should  be  taken  below  200  Hz . 

A7.  Source  depth  (in  m) 
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A8.  Maximum  range  of  track  (in  km) 

A9.  Parabolic  equation  flag  ;  0  means  that  the  standard  parabolic  equation 
Avill  be  used,  1  means  that  the  modified  parabolic  equation  will  be  used. 

AlO.  Spreading  flag  ;  0  means  that  the  source  is  a  point  source  in  a  3  D 
geometry  while  1  means  that  the  source  is  a  line  source  in  a  2-D  geometry  (i.e.  a 
horizontal  line  source  perpendicular  to  the  range-depth  plane). 

All.  Source  flag  ;  0  means  that  the  staner  will  be  the  Gaussian  source.  1  the 
starting  field  is  provided  by  the  user  in  a  file.  2  the  starter  is  the  Greene  source 
and  3  the  starter  is  the  Thomson  source. 

A12.  Half  beam  width  and  Beam  tilt  (up  is  positive)  with  respect  to  the  hor¬ 
izontal  (in  degrees)  for  source  flag  equal  to  0. 2.  3. 

A13.  Filename  containing  the  starting  field  when  it  is  provided  as  input 
(source  flag  equal  1). 

A 14.  Number  of  receiver  depths  and  number  of  selected  ranges  (each  one 
2  20)  where  transmission  loss  is  desired. 

A15.  Number  of  ranges  and  number  of  depths  for  the  contour  plot  (max¬ 
imum  dimension  of  100  for  each). 

A 16.  List  of  the  receiver  depth  array  (m)  in  increasing  order 

A17.  List  of  the  selected  ranges  (m)  in  increasing  order 

A 18.  Number  of  seaors  (the  environmental  data  must  be  repeated  for  each 
sector) 

Bl.  Water  depth  (in  m)  and  starting  range  of  the  sector  (in  km) 

B2.  Number  of  points  in  the  water  sound  velocity  profile 


B3.  List  of  pairs  depth,  sound  velocity  (in  m  and  nVs) 

B4.  Sediment  hei^i  in  m  and  density  (in  f/on*) 

B5.  Attenuation  in  the  sediments  ( in  ) 

B6.  Number  of  points  in  the  seditncnt  sound  velocity  profile 
B7.  List  of  the  couples  height,  sound  velocity  (in  m  and  m/s) 

B8.  Subboitom  density  (in  i/em  attenuauon  (in  dB  X) 

B9.  Subbottom  sound  velocity  (in  ) 

When  two  arguments  are  on  the  same  line,  be  sure  to  have  at  least  one 
space  (not  a  comma)  in  between.  The  couples  heighi/sound  velocity  must  be 
such  that  the  hrst  sound  velocity  value  corresponds  to  a  zero  height  and  the  last 
one  to  the  maximum  height  Also  note  that  heights  are  used  and  not  depths  for 
the  sediment  layer. 
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10.  Program  Output 

The  user  is  allowed  in  the  input  file  to  ask  for  five  different  outputs  : 


-  a  sector  information  file  that  summarizes  all  environmental  information 

-  a  transmission  loss  versus  range  file  that  contains  the  transmission  loss 
versus  range  for  each  receiver  depth  over  the  entire  range  of  track.  A  multi¬ 
channel  SIO  data  file  is  automatically  created  for  subsequent  processing  and 
plotting.  The  name  of  this  SIO  data  file  is  the  inputed  filename  with  the  "sio" 
extension.  The  first  channel  of  the  file  contains  the  ranges  and  the  following  ones 
(up  to  20)  correspond  to  a  receiver  depth  and  contain  the  transmission  loss. 

-  a  transmission  loss  versus  depth  that  contains  the  transmission  loss  versus 
depth  for  each  selected  range,  over  the  maximum  physical  depth.  A  multi- 
chaimel  SIO  data  file  is  automatically  created  for  subsequent  processing  and 
plotting.  The  name  of  this  SIO  file  is  the  inputed  filename  with  a  "sio"  extension. 
The  file  is  strucUired  in  pair  of  channels  where  the  2/7-1'*  and  the  2p"'  channel 
correspond  to  the  depth  and  the  transmission  loss  of  the  p"*  selected  range. 

-  a  SIO  data  file  to  feed  a  contour  routine.  This  file  has  a  number  of  chan¬ 
nels  equal  to  the  number  of  depth  samples  (A  15).  Each  channel  has  a  number  of 
points  equal  to  the  number  of  ranges  samples  (A  15).  The  file  is  an  equi-spaced 
grid  in  range  and  depth  of  transmission  loss  values. 

-  a  SIO  data  file  containing  the  pressure  field  where  the  first  channel  con¬ 
tains  the  real  part  of  the  pressure  and  the  second  one  the  imaginary  pan  of  the 
pressure.  This  output  requires  a  fixed  range  step  and  a  fixed  transform  size  N. 
Each  channel  is  made  of  sequential  2*^  point  segments  ;  Each  segment 
corresponds  to  a  range  sample. 


I 

I 
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11.  Examples 


11.1.  Isovelocity  Shallow  Water  (Range  Independent) 

The  first  example  is  a  shallow  water  run  taken  from  [1].  The  medium  is 
range  independent  with  isovelocity  water  and  bottom,  but  without  subbottom. 
The  source  used  is  the  standard  Gaussian  source  (half  bandwidth  of  35")  at  a 
depth  of  25  m .  The  frequency  is  50  Hi .  The  standard  parabolic  equation  is  used 
with  a  reference  sound  speed  equal  to  1500  »j/5.  No  phase  velocity  correction  is 
applied.  A  copy  of  the  script,  the  input  file  and  the  plot  control  files  are 
presented  below.  In  this  example  two  modes  are  propagating,  one  with  a 
stronger  attenuation  than  the  other.  We  find  the  same  transmission  loss  versus 
range  plot  as  in  [1].  This  mn  takes  15  minutes  on  a  Sun  Worksation  3/50  (Sun 
Workstation  is  a  registered  trademark  of  Sun  Microsystems,  Inc), 
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11^.  Wedge>Shaped  Ocean  with  Penetrable  Bottom 

The  second  example  tries  to  repeat  the  results  by  Jensen  and  Kuperman 
[17].  The  environment  consists  of  an  initial  5  km  stretch  of  200  m  depth  fol¬ 
lowed  by  a  bottom  slope  of  1.55*.  The  source  used  is  the  standard  Gaussian 
source  (the  half  width  is  not  given  in  [17]),  at  112  m  depth,  the  frequency  is  25  Hz . 
The  maximum  water  depth  is  200  m  with  a  sediment  layer  extending  to  a  depth  of 
800m. 

The  program  is  such  that  the  contour  plot  extends  in  range  from  the  source 
to  the  maximum  range,  and  in  depth  from  the  surface  to  the  maximum  depth. 
The  maximum  depth  is  the  maximum  sum  of  the  water  depth,  sediment/bottom 
and  subbottom  heights. 

The  transmission  loss  versus  depth  extends  in  depth  from  the  surface  to  the 
maximum  depth  like  in  the  contour  plot,  and  all  mesh  points  axe  used.  The 
transmission  loss  versus  depth  is  interpolated  in  range  when  the  selected  range  is 
not  a  range  mesh  point. 

The  contour  plot  is  presented  below.  One  finds  good  agreement  with  [17] 
although  the  contour  routines  are  different.  A  copy  of  the  input  file  as  well  as  the 
script  and  plot  control  files  and  the  sector  information  output  file  are  following. 
This  job  takes  50  minutes  to  run  on  a  Sun  Workstation  3/50  (Sun  Workstation  is  a 
registered  trademark  of  Sun  Microsystems,  Inc). 
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lU.  Deep  Water  Example  I  (Range  Independent) 

This  example  is  of  deep  water  propagation  over  a  range  of  60  km  ,the  water 
depth  is  2000  m ,  the  frequency  is  100  Hz  and  the  source  is  a  T  narrow  and  8"  tilted 
Gaussian  source  at  250  m  depth  [1].  One  finds  once  again  agreement  between  the 


MPL’s  PE  program  and  [1].  Script,  input,  plot  control  files  and  plots  follow. 
This  run  takes  30  minutes  on  a  Sun  Workstation  3/50  (Sun  Workstation  is  a 
registered  trademark  of  Sun  Microsystems,  Inc). 
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GUOBAL  PLOT  COWTROL  PARAMETERS  0,  numsyn-  symbol  number  for  point  plotting  (37) 
nun^lt'  number  of  data  files  plotted  on  one  plot  (2)  0.15,  symsiz-  symbol  size  in  inches  (38) 
xlen-  X  axis  length  in  inches  (3)  0.0,  xoff-  x  offset  in  user  units  (39) 
noaxe  -  axis  suppression.  (0)  supress,  (1)  box,  (2)  lower,  (3)  mid  (4)  0.0,  yoff-  y  offset  in  user  units  (40) 
nsf-  san^ling  frequency  flag  (S)  0,  ndash  -  number  of  dashed  lines  for  this  plot  (41) 
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11.4.  Deep  Water  Example  n  (Range  Independent) 

This  last  example  is  of  deep  water  propagation  at  100  Hz  in  a  range  indepen¬ 
dent  medium  [1].  The  source,  250  meter  deep  is  a  Gaussian  source  with  a  T 
width  and  a  8®  tilt.  The  maximum  water  depth  is  2000  m  with  a  100  m  sediment 
layer  and  a  subbottom.  The  script  produces  a  plot  of  transmission  loss  versus 
range,  a  plot  of  transmission  loss  versus  depth  and  a  partial  contour  plot.  One  can 
see  agreement  with  [1].  Script,  input,  plot  control  files  and  plots  follow.  This  run 
takes  30  minutes  on  a  Sun  Worstation  3/50  (Sun  Workstation  is  a  registered 
trademark  of  Sun  Microsystems,  Inc). 
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Example  4  :  Deep  Water  (0,0,100,0)  **  TILTED  BEAM 
Freq'T:  100  Hz,  ^SD:  250  m,  RR:  40  l<ni,  GAUSS  (2,8) 
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